library(here)
library(readxl)
#library(papeR)
#library(outliers)
library(kableExtra)
#library(DataExplorer)
library(lubridate)
library(forecast)
library(nlme)
#library(nortest)
library(ggfortify)
library(dygraphs)
#library(seasonal)
#library(seasonalview)
library(nonlinearTseries)
library(fNonlinear)
library(fGarch)
library(TSA)
library(tsDyn)
library(tidyverse)
Rutas
raw_data <- here("data", "raw")
interim_data <- here("data", "interim")
final_data <- here("data", "processed")
Leyendo base de datos
base <- read_excel(paste0(raw_data,"/Base Datos.xlsx"),
col_types = c("text", "numeric", "numeric"))
head(base)
Extrayendo la base de dolares y de dolares
dolares <- data_frame(date = base$`Activo neto`,
dol = as.double(base$USD)) %>%
mutate(date = ymd(paste0(date, "-01"))) %>%
mutate(year = as.factor(year(date)),
month = as.factor(month(date)))
head(dolares)
Convirtiendo los datos a series de tiempo
# dolares
# 6 periodos para utilizar a modo de validación
dolares_val <- dolares %>%
slice_tail(n = 6)
# Serie completa para el analisis exploratorio
dolares_full <- dolares
# 10 años de datos para modelar
dolares <- dolares %>%
filter(date>"2007-12-01" & date<="2020-01-01")
# Objetos ts tanto para la serie completa como para la serie de modelado
dolares_ts_full <- ts(dolares_full$dol, start = c(2001,2), frequency = 12)
dolares_ts <- ts(dolares$dol, start = c(2008,1), frequency = 12)
dolares_ts_val <- tail(dolares_ts_full, 6)
# Null hypothesis: Linearity in "mean"
tnnTest(dolares_ts, lag = 1, title = NULL, description = NULL)
Title:
Teraesvirta Neural Network Test
Test Results:
PARAMETER:
lag: 1
m|df: 2
t-lag-m|df: 142
STATISTIC:
Chi-squared: 6.55
F: 3.2808
P VALUE:
Chi-squared: 0.03782
F: 0.04047
Description:
Thu Sep 10 22:21:31 2020 by user:
La hipótesis nula del Teraesvirta Neural Network Test es que la media de la serie es lineal. Al rechazarse la hipótesis con una confianza del 95% se puede decir que la serie es no lineal.
options(max.print=1000000)
rqa.analysis=rqa(time.series = dolares_ts, embedding.dim=1, time.lag=1,radius=0.01,lmin=2,vmin=2,do.plot=TRUE,distanceToBorder=2)
# Una forma visual de empezar a revisar si existen o no clusters de volatilidad
dolares_ts_nd <-diff(log(dolares_ts))
dolares_ts_nd2<-dolares_ts_nd-mean(dolares_ts_nd) # Es el cambio relativo ajustado por la media en el tipo de cambio
#dolares_ts_nd2
dolares_ts_nd3<-dolares_ts_nd2^2 # Medida de la volatilidad. Al ser una cantidad al cuadrado, su valor ser? alto en periodos en que se experimenten grandes cambios y comparativamente peque?o cuando sucedan cambios modestos en los precios de dichos bienes.
plot(dolares_ts_nd3)
En este gráfico de la serie ajustada por la media y elevada al cuadrado se pretende observar la volatilidad de la misma. Al ser una cantidad al cuadrado, cuando su valor ses alto en indica que se experimenten grandes cambios y comparativamente pequeño cuando sucedan cambios modestos en los precios de dichos bienes. Es así que es claro que después del 2008, alrededor del 2010 y en el 2013 es donde se presentan los mayores cambios comparativos y es efectivamente en donde se identifican crisis económicas que llevaron a la gente a buscar estas opciones de inversión.
plot(dolares_ts_nd,type="l"); abline(h=0)
qqnorm(dolares_ts_nd); qqline(dolares_ts_nd)
acf(as.vector(dolares_ts_nd))
pacf(as.vector(dolares_ts_nd))
#Graficos para corroborar independencia(ruido blanco) que es diferente de correlaci?n (medida de dependencia lineal).
#Ho:residuos son independientes
acf(dolares_ts_nd^2)
pacf(dolares_ts_nd^2)
acf(abs(dolares_ts_nd))
pacf(abs(dolares_ts_nd))
En este caso algunas estacas se salen (algunas autocorrelaciones son significativas) y por tanto los rendimientos no son independientes ni identicamente distribuidos. Las autocorrelaciones significativas de los rendimientos al cuadrado o en términos absolutos reflejan la existencia de agrupamiento de volatilidad.
#McLeod.Li (Box-Ljung) test muestra una evidencia fuerte de heterocedasticidad condicional(p-value significativo).
TSA::McLeod.Li.test(y=dolares_ts_nd)
Además, a partir del test de McLeod.Li (Box-Ljung) se muestra evidencia fuerte de heterocedasticidad condicional ya que varios p-value son significativos.
garch_mod <- garchFit(~arma(0,0)+garch(1,1), data=dolares_ts_nd,include.mean = FALSE)
Series Initialization:
ARMA Model: arma
Formula Mean: ~ arma(0, 0)
GARCH Model: garch
Formula Variance: ~ garch(1, 1)
ARMA Order: 0 0
Max ARMA Order: 0
GARCH Order: 1 1
Max GARCH Order: 1
Maximum Order: 1
Conditional Dist: norm
h.start: 2
llh.start: 1
Length of Series: 144
Recursion Init: mci
Series Scale: 0.07852679
Parameter Initialization:
Initial Parameters: $params
Limits of Transformations: $U, $V
Which Parameters are Fixed? $includes
Parameter Matrix:
Index List of Parameters to be Optimized:
omega alpha1 beta1
2 3 5
Persistence: 0.9
--- START OF TRACE ---
Selected Algorithm: nlminb
R coded nlminb Solver:
0: 205.61475: 0.100000 0.100000 0.800000
1: 205.43217: 0.109958 0.0990399 0.807857
2: 205.21550: 0.109406 0.0867192 0.804739
3: 204.82498: 0.125244 0.0698920 0.815383
4: 204.72240: 0.135542 0.0209029 0.824502
5: 204.48753: 0.184347 0.00695293 0.820955
6: 204.33288: 0.211453 0.0115362 0.778136
7: 204.29614: 0.172368 0.0155603 0.810467
8: 204.28448: 0.156327 0.0208586 0.821820
9: 204.27178: 0.122608 0.0206070 0.855240
10: 204.27069: 0.122959 0.0207417 0.855573
11: 204.27011: 0.123074 0.0202577 0.855642
12: 204.26812: 0.125112 0.0185096 0.855599
13: 204.26796: 0.124847 0.0177023 0.856713
14: 204.26795: 0.125316 0.0177612 0.856194
15: 204.26795: 0.125235 0.0177596 0.856275
16: 204.26795: 0.125236 0.0177593 0.856275
Final Estimate of the Negative LLH:
LLH: -162.1135 norm LLH: -1.125788
omega alpha1 beta1
0.0007722607 0.0177593449 0.8562746624
R-optimhess Difference Approximated Hessian Matrix:
omega alpha1 beta1
omega -87090880.9 -520359.912 -534764.923
alpha1 -520359.9 -3752.424 -3245.468
beta1 -534764.9 -3245.468 -3319.657
attr(,"time")
Time difference of 0.01039481 secs
--- END OF TRACE ---
Time to Estimate Parameters:
Time difference of 0.1755028 secs
summary(garch_mod)
Title:
GARCH Modelling
Call:
garchFit(formula = ~arma(0, 0) + garch(1, 1), data = dolares_ts_nd,
include.mean = FALSE)
Mean and Variance Equation:
data ~ arma(0, 0) + garch(1, 1)
<environment: 0x7ff522bb3a80>
[data = dolares_ts_nd]
Conditional Distribution:
norm
Coefficient(s):
omega alpha1 beta1
0.00077226 0.01775934 0.85627466
Std. Errors:
based on Hessian
Error Analysis:
Estimate Std. Error t value Pr(>|t|)
omega 0.0007723 0.0010342 0.747 0.455
alpha1 0.0177593 0.0417714 0.425 0.671
beta1 0.8562747 0.1764943 4.852 1.22e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Log Likelihood:
162.1135 normalized: 1.125788
Description:
Thu Sep 10 22:23:35 2020 by user:
Standardised Residuals Tests:
Statistic p-Value
Jarque-Bera Test R Chi^2 2.496111 0.2870624
Shapiro-Wilk Test R W 0.9929834 0.7046679
Ljung-Box Test R Q(10) 13.40665 0.2018155
Ljung-Box Test R Q(15) 21.8283 0.1123856
Ljung-Box Test R Q(20) 31.40617 0.05005161
Ljung-Box Test R^2 Q(10) 6.433964 0.7775818
Ljung-Box Test R^2 Q(15) 10.35613 0.7967609
Ljung-Box Test R^2 Q(20) 12.30276 0.9052611
LM Arch Test R TR^2 6.194549 0.9059588
Information Criterion Statistics:
AIC BIC SIC HQIC
-2.209909 -2.148038 -2.210754 -2.184768
acf(residuals(garch_mod)^2)
pacf(residuals(garch_mod)^2)
#1. Test de Portmanteu para residuos estandarizados al cuadrado donde la Ho es que los residuos no est?n correlacionados.
try(
gBox(garch_mod,method="absolut", plot = T)
)
Error in model$call : $ operator not defined for this S4 class
A partir de esta prueba y sumado a la inspección visual se determina que los residuos del modelo Garch estan correlacionados pues el p-value es mayor al punto de corte (0.05) aunque las estacas en los gráficos no son significativas.
fGarch::predict(garch_mod, n.ahead = 6,mse="uncond",plot=TRUE, crit_val=2)
NA
Por último, se observa que el ajuste del modelo en términos de predicción es muy deficiente.
nn_mod<-nnetar(dolares_ts_nd)
summary(nn_mod)
Length Class Mode
x 144 ts numeric
m 1 -none- numeric
p 1 -none- numeric
P 1 -none- numeric
scalex 2 -none- list
size 1 -none- numeric
subset 144 -none- numeric
model 20 nnetarmodels list
nnetargs 0 -none- list
fitted 144 ts numeric
residuals 144 ts numeric
lags 2 -none- numeric
series 1 -none- character
method 1 -none- character
call 2 -none- call
acf(residuals(nn_mod)[!is.na(residuals(nn_mod))]^2)
pacf(residuals(nn_mod)[!is.na(residuals(nn_mod))]^2)
#1. Test de Portmanteu para residuos estandarizados al cuadrado donde la Ho es que los residuos no est?n correlacionados.
gBox(nn_mod,method="absolut", plot = T)
A partir de esta prueba y sumado a la inspección visual se determina que los residuos del modelo Garch estan correlacionados pues el p-value es mayor al punto de corte (0.05) y un par de estacas en los gráficos son significativas.
pred_nn_mod<-forecast::forecast(nn_mod,level = c(95), h=6, bootstrap=TRUE, npaths=10000)
pred_nn_mod
plot(pred_nn_mod)
Finalmente, se observa que aunque mejor que el modelo Garch, el modelo de redes neuronales tampoco tienen un buen ajuste en las predicciones.
dimension = estimateEmbeddingDim(dolares_ts_nd, time.lag=1, max.embedding.dim=15,threshold=0.95, do.plot=TRUE)
aar_mod <- aar(dolares_ts_nd, m=dimension)
summary(aar_mod)
Non linear autoregressive model
AAR model
Family: gaussian
Link function: identity
Formula:
y ~ s(V1.0, bs = "cr") + s(V1..1, bs = "cr") + s(V1..2, bs = "cr") +
s(V1..3, bs = "cr") + s(V1..4, bs = "cr") + s(V1..5, bs = "cr") +
s(V1..6, bs = "cr") + s(V1..7, bs = "cr") + s(V1..8, bs = "cr")
Estimated degrees of freedom:
8.29 1.00 1.94 1.00 2.95 1.00 2.26
1.00 1.00 total = 21.44
GCV score: 0.006481423
Residuals:
Min 1Q Median 3Q Max
-0.16838301 -0.04239723 0.00050029 0.04749853 0.18649028
Fit:
residuals variance = 0.0043, AIC = -621, MAPE = 168.6%
Family: gaussian
Link function: identity
Formula:
y ~ s(V1.0, bs = "cr") + s(V1..1, bs = "cr") + s(V1..2, bs = "cr") +
s(V1..3, bs = "cr") + s(V1..4, bs = "cr") + s(V1..5, bs = "cr") +
s(V1..6, bs = "cr") + s(V1..7, bs = "cr") + s(V1..8, bs = "cr")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0064205 0.0063550 1.0103 0.3145
Approximate significance of smooth terms:
edf Ref.df F p-value
s(V1.0) 8.2911 8.7922 1.1994 0.276636
s(V1..1) 1.0000 1.0000 0.9257 0.338013
s(V1..2) 1.9394 2.4627 0.7609 0.413436
s(V1..3) 1.0000 1.0000 0.0034 0.953689
s(V1..4) 2.9511 3.6824 1.5087 0.189160
s(V1..5) 1.0000 1.0000 0.4893 0.485682
s(V1..6) 2.2586 2.8659 4.3585 0.006241 **
s(V1..7) 1.0000 1.0000 0.0313 0.859901
s(V1..8) 1.0000 1.0000 0.2346 0.629058
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.123 Deviance explained = 25.7%
GCV = 0.0064814 Scale est. = 0.0054521 n = 135
plot(aar_mod)
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
e_aar_mod <- residuals(aar_mod)
plot(e_aar_mod)
e_aar_mod <- e_aar_mod[!is.na(e_aar_mod)]
acf(e_aar_mod)
pacf(e_aar_mod)
AIC(aar_mod)
[1] -620.6911
mse(aar_mod)
[1] 0.004299555
MAPE(aar_mod)
[1] 1.686104
#fitted(aar_mod)
#coef(aar_mod)
pred_aar <- predict(aar_mod, n.ahead=6)
autoplot(ts(c(dolares_ts_nd, pred_aar), start = start(dolares_ts_nd), frequency = frequency(dolares_ts_nd)) )
Se observa que las medidas de ajuste (MSE y MAPE) no son malas y la predicción para los 6 periodos es mejor que el modelo Garch pero hay que comparar con los demás para determinar cual obtiene mejor ajuste.
star_mod <- star(dolares_ts_nd, mTh=c(0,1), control=list(maxit=10000))
Testing linearity... p-Value = 0.1801151
The series is linear. Using linear model instead.
summary(star_mod)
Non linear autoregressive model
AR model
Coefficients:
const phi.1 phi.2
0.007107589 -0.017957086 -0.049583561
Residuals:
Min 1Q Median 3Q Max
-0.196202131 -0.050501541 -0.005116169 0.049886545 0.283224845
Fit:
residuals variance = 0.005981, AIC = -731, MAPE = 125.4%
Coefficient(s):
Estimate Std. Error t value Pr(>|t|)
const 0.0071076 0.0066029 1.0764 0.2836
phi.1 -0.0179571 0.0838235 -0.2142 0.8307
phi.2 -0.0495836 0.0835235 -0.5936 0.5537
plot(star_mod)
NA
NA
e_star_mod <- residuals(star_mod)
plot(e_star_mod)
e_star_mod <- e_star_mod[!is.na(e_star_mod)]
acf(e_star_mod)
pacf(e_star_mod)
AIC(star_mod)
[1] -731.1572
mse(star_mod)
[1] 0.005981122
MAPE(star_mod)
[1] 1.253908
pred_star <- predict(star_mod, n.ahead=6)
autoplot(ts(c(dolares_ts_nd, pred_star), start = start(dolares_ts_nd), frequency = frequency(dolares_ts_nd)) )
Los indicadores del modelo Star son mejores que los del modelo Aar, tanto el AIC, como el MSE y MAPE. Los resultados del pronóstico parecen ser peores que el modelo AAR y es por esta razón que se procede a revisar los pronósticos de todos los modelos no lineales para determinar el mejor.
Una vez teniendo los modelos de suavizamiento exponencial, de regresión y Box-Jenkings (ARIMA) se procede a compararlos en términos de ajuste visual y de indicadores de ajuste para determinar cual es el mejor modelo que pronostique el número de muertes por accidentes de tránsito en Costa Rica.
# Prediccion redes neuronales
f_mod1 <-
exp(
log(dolares_ts[145]) + cumsum(pred_nn_mod$mean)
)
f_mod1_ts = ts(
f_mod1,
frequency = 12,
start = c(2020, 2),
end = c(2020, 7)
)
# Prediccion modelo aar
f_mod2 <-
exp(
log(dolares_ts[145]) + cumsum(pred_aar)
)
f_mod2_ts = ts(
f_mod2,
frequency = 12,
start = c(2020, 2),
end = c(2020, 7)
)
# Prediccion modelo Star
f_mod3 <-
exp(
log(dolares_ts[145]) + cumsum(pred_star)
)
f_mod3_ts = ts(
f_mod3,
frequency = 12,
start = c(2020, 2),
end = c(2020, 7)
)
todos_preds <- cbind(
Serie = dolares_ts,
Real = dolares_ts_val,
Prediccion1 = f_mod1_ts,
Prediccion2 = f_mod2_ts,
Prediccion3 = f_mod3_ts
)
# Graficamos
dygraph(todos_preds, main = "Predicción todos los modelos") %>%
dySeries("Serie", label = "Cantidad") %>%
dySeries("Real", label = "Real") %>%
#dySeries("Prediccion1", label = "SES") %>%
#dySeries("Prediccion2", label = "Holt") %>%
dySeries("Prediccion1", label = "Red_neuronal") %>%
dySeries("Prediccion2", label = "AAR") %>%
dySeries("Prediccion3", label = "STAR") %>%
dyAxis("x", label = "Meses") %>%
dyAxis("y", label = "Accidentes") %>%
dyOptions(colors = RColorBrewer::brewer.pal(7, "Set1")) %>%
dyRangeSelector()
acc_todos <- tibble(
Metodo = c("Red_neuronal", "AAR", "STAR"),
RMSE = round(
c(
forecast::accuracy(f_mod1_ts, dolares_ts_val)[2],
forecast::accuracy(f_mod2_ts, dolares_ts_val)[2],
forecast::accuracy(f_mod3_ts, dolares_ts_val)[2]
),
3
),
MAE = round(
c(
forecast::accuracy(f_mod1_ts, dolares_ts_val)[3],
forecast::accuracy(f_mod2_ts, dolares_ts_val)[3],
forecast::accuracy(f_mod3_ts, dolares_ts_val)[3]
),
3
),
MAPE = round(
c(
forecast::accuracy(f_mod1_ts, dolares_ts_val)[5],
forecast::accuracy(f_mod2_ts, dolares_ts_val)[5],
forecast::accuracy(f_mod3_ts, dolares_ts_val)[5]
),
3
)
)
acc_todos %>%
mutate_if(is.numeric, function(x) {
cell_spec(x, bold = T,
color = spec_color(x, end = 0.9, direction = -1),
font_size = spec_font_size(x, begin = 12,end = 14, scale_from = NA))
}) %>%
kable(escape = F, caption = "Medidas de ajuste (Todos los modelos)", digits = 2) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
| Metodo | RMSE | MAE | MAPE |
|---|---|---|---|
| Red_neuronal | 174.818 | 141.332 | 10.069 |
| AAR | 146.611 | 119.596 | 8.503 |
| STAR | 169.26 | 138.179 | 9.855 |
NA
Finalmente, se obtienen los valores de pronostico de 3 modelos a contrastar: Red neuronal, AAR y STAR. Se puede apreciar que el de red neuronal fue el mas alejado de los valores reales. Entre los otros modelos se puede observar que tienen, en términos gráficos, resultados similares aunque difieren al final del periodo de prueba donde el AAR ajusta mejor. Moviendonos al análisis de los indicadores de ajuste se puede ver que el ARR fue el modelo que obtuvo los mejores indicadores y es entonces que se selecciona como el mejor modelo no lineal para predecir el número de autos importados por mes en Costa Rica.